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ABSTRACT 


The near-field, explosive dispersal of a liquid into air has been 
explored using a combination of analytical and numerical models. The 
near-field flow regime is transient, existing only as long as the 
explosive forces produced by the detonation of the burster charge 
dominate or are approximately equal in magnitude to the aerodynamic drag 
forces on the liquid: The near-field model provides reasonable initial 
conditions for the far-field model, which is described in a separate 
report. The near-field model consists of the CTH hydrodynamics code and 
a film instability model. In particular, the CTH hydrodynamics code is 
used to provide initial temperature, pressure, and velocity fields, and 
bulk material distribution for the far-field model. The film 
instability model is a linear stability model .for a radially expanding 
fluid film, and is used to provide a lower bound on the breakup time and 
an upper and lower bound on the initial average drop diameter for the 
liquid following breakup. Predictions of the liquid breakup time and 
the initial arithmetic average drop diameter from the model compare 


favorably with the sparse experimental data. 
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1. INTRODUCTION 


. A fuel-air explosive device mixes a liquid or solid fuel with air to 


produce a detonable mixture. The archetypical fuel-air explosive device 
is cylindrical, with a central, cylindrical burster charge surrounded by 


a cylindrical annulus of liquid fuel. When the central burster charge 


is detonated, the casing is fragmented, and the fuel is propelled 


explosively outward and breaks into drops. Further expansion and 
turbulent mixing create an aerosol cloud that is roughly toroidal in 


shape. Detonating the resulting cloud with a second high explosive — 


charge can produce a blast wave which can inflict considerable damage on 
blast- sensitive targets ‘such as buildings and vehicles. 


Current fuel-air explosive sevice use liquid fuels with high vapor 


pressures. Such fuels pose a,/significant safety hazard, as leaking fuel 
can evaporate and form a detonable fuel-air ‘mixture in a confined space 
(such as on board,a ship). Such devices thus require hypergolic 
storage, which limits their availability and utility. 


Low vapor pressure liquid fuels and solid fuels, and slurried and gelled 
fuels promise greater patety, if they can be Peraeoly setoneted after 


dispersal. : : | oo 


This includes understanding the mechanism(s) of fuel dispersal and 
consequently which factors influence the final cloud shape and fuel 


- Improving the reliability of detonation of a fuel-air cloud depends on 


understanding the mechanisms important in the detonation of aerosols. 


distribution produced by dispersal from a. given device. 


herb iwst: yeeta is a function of Pieieee reaction chemistry, which 
depends on the physical state of the fuel including particle size and 
concentration. The latter are determined by the dispersal process. 

Thus the outcome of the dispersal process impacts the reliability of 
detonation of the fuel-air cioud, and also the blast yield. 
Consequently, an understanding |of the process of explosive dispersal is 
important in designing improved fuel-wir explosive couicen: . 


The Enhanced Blast Munitions Program has addressed how to provide the 
technology base, that would enable the design and development of safe, 
reliable, high-yield, and high-efficiency blast weapons based on fuel- 
air explosives. Within that program, the goal of the dispersal modeling 
work Was to create a model to simulate the formation of a fuel-air 
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aerosol cloud, from the initial device to the point of detonation of the 
fuel-air cloud. Ultimately, the model was intended to include liquid 
and solid fuels, and devices. of various sizes and shapes. 


The process of explosive dispersal can be divided into three regimes, 


based on the relative magnitudes of the explosive and aerodynamic forces 


acting on the fuel [1]. In the ejection regime, explosive forces from 
the detonation of the burster charge dominate the aerodynamic forces. 
This regime is succeeded. by the transition regime, in which the 
explosive and aerodynamic forces are of approximately equal magnitude. 

In both the ejection and transition regimes, the fuel concentration in 
the evolving cloud is relatively high and strong shock waves reverberate 
through the cloud. The transition regime is succeeded by the expansion 
regime, in which aerodynamic forces dominate the explosive forces, the 
fuel concentration in the cloud is relatively low, and the shock waves 
have decayed to insignificance. 


Previous dispersal modeling efforts for fuel-air explosives have. 
included both linear-regressions of vast amounts of experimental data 
(e.g., [1, 2, 3]), and numerical solutions to the Navier-Stokes 
equations (e.g., [3, 4, 5]). However, none of these previous efforts 
modeled the entire dispersal process, and in particular no attempt was 
made to model the initial breakup of liquids into drops [6]. 


For the purposes of our dispersal modeling effort, the ejection and | 
transition regimes described above were combined into a regime called 

the near-field regime. This definition reflects the fact that less is 
known about the physical processes occurring in the ejection and 
transition regimes than in the expansion regime. To complement the new 


_ terminology, the expansion regime was renamed the far-field regime. 


The general approach employed in our dispersal modeling effort in each 
of the near-field and far-field regimes was to 


* use existing numerical tools where possible, 
* consider solid and liquid fuels separately, and ™ 
* endeavor to understand the underlying physics. 


Consonant with this general approach, and early in the dispersal 
modeling effort, the KIVA hydrodynamics code [7-10] was selected as the 


W misqua t ryt mI 1 1h 


far-field modeling tool. KIVA is a one-, two-, or three-dimensional. 
finite difference code for modeling transient fluid flows with 
chemically-reacting fuel sprays. The fuel drops are modeled 
stoichastically, using parcels to represent a number of physical 
particles with identical properties. The desired initial drop size 
distribution can be chosen by the user. KIVA incorporates models for 
the evaporation, aerodynamic breakup, collision and coalescence, and 
turbulent motion of liquid fuel drops. It also incorporates both a k-e 
turbulence model and a subgrid turbulence model. (Although turbulence 
models are not well-developed for compressible flows in unconfined 
geometries, the presence of these models in KIVA was considered 
sufficient for the initial phases of the development of the dispersal 
models.) Hence it was considered to be a good choice for simulating the 
formation of aerosol clouds. _ | 


However, KIVA is inaccurate for large fuel concentrations, and is unable 
to model well the strong shock waves present during the near-field 
regime of the dispersal process. 


In the context of our dispersal modeling work, the near-field regime 
extends in time from the detonation of the burster charge in the fuel- 
air device to a time when the average fuel concentration in the 
developing aerosol cloud is small enough and the shock waves from the 
burster charge are weak enough for KIVA to be a valid model. The near- 
field regime thus includes the processes of 


* fragmentation of the outer canister enclosing the fuel, 
e the early expansion of the fuel as a coherent shell, and 
e the breakup of the shell into a field of drops. 


The spatial extent of the near-field regime depends on the device 
dimensions. A conservative estimate of the maximum radial extent of the 
near-field regime, based on concentration only, is shown in Figure 1, as 
a function of the outer radius of the device. The radial extent of the 
near-field regime is over-estimated in Figure 1. This figure shows the 
estimated outer radius of the fuel air cloud at the time when the near- 
field regime vanishes. 


This report describes the near-field dispersal modeling effort. The 
far-field dispersal modeling effort is described in [23]. The coupled 
near-field, far-field model is described in [22]. 
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Figure 1: The maximum radial extent of the near-field regime increases 
with increasing dimensions of a fuel-air explosive device 
(initial device inner radius = 1 cm). 
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The goal of the wear-field modeling effort is to provide reasonable 
{initial conditions for the far-field model. This includes specifying 
appropriate initial 


° temperature fields, 
: pressure fields, 
¢ velocity fields, 
— @ bulk naterial distributions, 
© average particle sizes, and. 
e particle size dveueipueicas? 


Modeling of formation of the fuel-air aerosol cloud is hindered by the 
sparsity of reliable experimental data. Most of the data reported in 
the literature are the average rate of cloud growth (radius and height, 
or volume) for various fuels and sizes and geometries of devices. 
Various attempts have been made to measure particle sizes in fuel-air 
clouds, but with mixed success [1, 3, 11, 12, 13, 14]. Almost all the 
data available are for the fully developed cloud (i.e., in the far-field 
regime, as defined here). The only data available for the near-field 
regime were reported by Samirant et al. [12]. 


The numerical model used to predict the temperature, pressure, and 
velocity fields, and the bulk material distribution, will be described 
first. Next the complementary analytical model will be described; this 
model was developed to predict the liquid breakup time and the 
arithmetic average drop diameter at breakup. The resulting combined 
near-field model will be illustrated, and recommendations for further 
work will be discussed. Some details of the analytical model are 
presented in Appendix 4, and two energy-based breakup models are 
discussed in Appendices 2 and 3. 
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2. NUMERICAL MODELS 


Numerical modeling of the breakup of a mass of liquid into a field of 
drops is intrinsically difficult, owing to 


* the many length scales involved, 
* the important role played by interfacial physics, 
* the need to ge une: location of material phases, and 


* the requirement that the physics of Liquids and their 
vapors be treated simultaneously. 


These factors place stringent limitations on any numerical modeiing 
tools used to simulate the explosive breakup of a liquid. 


The presence of multiple length scales complicates any numerical 
simulation because all the significant length scales must be adequately 
resolved. In a finite difference code, very fine computational grids 
are required. This greatly reduces the maximum allowed time step for 


any explicit solution scheme, and hence greatly increases computation 


costs. 


When the bulk liquid breaks up into a field of drops, there is an 
enormous increase in the ratio of surface area to volume for the liquid. 
Interfacial effects can no longer be ignored, i.e., surface tension 
becomes important. Certainly the representation of the surfaces becomes 
very important. This places several constraints on any numerical 
simulation tool that is to be used to simulate liquid breakup: 


* there must be little numerical diffusion, so that 
simulated interfaces remain distinct, 


¢ there must be an accurate representation of the interface, 
and 


* surface tension must be modeled correctly. 
The requirement of little numerical diffusion restricts the types of 


solution algorithms that can be used. The requirement for a good 
representation of the interface requires that some appropriate interface 


tracker be used. Surface tension is not modeled well in most hydro- 


dynamics codes. 


As the shock waves from the detonation of the central high explosive 
burster charge reverberate in the liquid mass, some of the liquid may 
vaporize [13]. Hence regions of a liquid and its vapor must be tracked 
independently. This is difficult to do, especially on the fine scales 
of cavitation seen in explosively dispersed liquids, and is not done in 
most hydrodynamics codes. 


In addition, the coexistence. of regions of condensed liquid and vapor 
presents problems for many hydrodynamics codes. Codes that work well 
for incompressible fluids may not work well at all for eompreseivre 
ones, and vice versa. 


Consonant with the strategy of using existing computational tools where 
possible, a survey of existing numerical tools was conducted. The 
survey included two- and three-dimensional codes, Lagrangian, Eulerian, 
and arbitrary Lagrangian-Eulerian (ALE) codes, vortex codes, particle- 
in-cell (PIC) codes, free~Lagrangian codes, finite difference codes, and 
finite element codes. Codes and algorithms were evaluated on the basis 
of the criteria discussed above. 


Following this survey, the CTH [15] hydrodynamics code was selected be- 
cause | 


* it was readily available, 


e there was little evidence of numerical diffusion in test 
‘Simulations of the dispersal of a liquid into air, 


e an explicit interface tracker is employed, 
* sophisticated equations of state can be employed, and 
e detonations of high explosives can be simulated. 


The CTH hydrodynamics code is an Eulerian finite difference code which 
uses a van Leer monotonic advection scheme that is formally second-order 
accurate. It can model multi-dimensional, multi-material shcck-wave 
physics. Simulations can be performed in one-, two-, or three-spatial 
dimensions, thus allowing the simulation of fuel-air explosive devices 
of various shapes. Models are included for treating the strength, 
fracture, and distension of various materials, and the detonation of 
high-explosives. 


I 


CTH leeks any model for surface tension; however, this was Prue for most 
codes ences in the survey. 


A typical dispereat simulation is shown in Figures 2-7 for liquid decane 
being dispersed into air by a cylindrical device. The central burster 
charge is the high explosive PBX9010 and the fuel mass to burster mass 
ratio was 100:1. The initial dimensions of one quarter of the model 
device are given in Figure 2. The CTH input files used to generate this 
simulation are presented in Appendix 1. This simulation models the 
liquid-fuel jug tests conducted at the Naval Weapons Center, China Lake, 
CA. Decane was chosen to represent transportation fuels, some of which 
were considered as fuels for improved fvel-air explosives. 


In each of Figures 2-7 only one quarter of the symmetrical problem 
domain is shown. The full domain is obtained by reflecting the region 
shown in both the axial and radial axes. 


Figure 3 shows the evolution of the material interfaces from 0 to 2 
msec. The decane expands radially and axially as a more or less 
coherent film. (The axial expansion could be inhibited by placing a 
metal plate on top of the device, instead of using a metal plug. This 
is often done in jug tests at the Naval Weapons Center, China Lake, CA.) 
The isopycnics (constant density contours) at 2 msec are shown in Figure 
4. Although they suggest that the liquid decane is not a coherent film, 
this is a computational artifact resulting from the finite size of the 
computational cells. 


In Figure 5, the steep thermal gradient across the decane film at 2 msec 
is evident. The temperature changes from approximately 1800 K at the 
inner surface of the film (adjacent to the detonation product gases) to 
300 K at the outer euneace of the film, over a distance on the order of 
t ce 


Figure 6 indicates that the pressure of the detonation product gases is 
less than ambient at 2 msec; the cloud continues to expand outward due 
to inertia. This is seen by the velocity vectors at 2 msec in Figure 7 
in which incipient vortices are visible. 


CTH is suitable for simulating the temperature, pressure, and velocity 
fields, and the bulk material distribution in the explosive dispersal of 
a liquid. However, owing to the absence of a surface tension model, it 
is unable to provide much useful information regarding the breakup of 
the liquid into drops: the time of breakup or the arithmetic average 
drop diameter immediately after breakup. Analytical models were 
investigated to remedy this deficiency. | 
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8. ANALYTICAL MODELS 


Owing to the difficulties discussed above, existing numerical models are 
unable to provide useful estimates of the initial arithmetic average 
drop diameter in the explosive breakup of a mass of liquid. Thus 
several analytical models were explored to provide the necessary 
estimates. The analytical models will provide some useful bounds on the 
drop diameters. | | 


Three analytical models were explored. The first model is a simple one 
based on minimizing surface energy. It provides estimates of either the 
breakup time (if the initial drop diameter is known) or the initial drop 
diameter (if the breakup time is known). This model is described in 
Appendix 2. | 


The second model is one devised by D. E. Grady [16-19]. It is also 
based on energy concepts. It provides estimates of both the breakup 
time and the initial arithmetic average drop diameter, and is described 
in Appendix 3, However, the model is based on the concept of liquid 
spall, which has not been experimentally observed [1, 20], and its 
predictions do not compare well with the (sparse) experimental data. 


The third analytical model is based on the linear instability of a 
radially expanding fluid film. It provides estimates of the breakup 
time and the initial arithmetic average drop diameter. Owing to the 
assumptions made in constructing the model, it provides a lower bound on 
the breakup time, and both upper and lower bounds on the initial 
arithmetic average drop diameter. Its predictions compare favorably 
with the (sparse) experimental data, and hence it is the model 
recommended by this study. 


The film instability model is described below. 


3.1 Description of the Film Instability Model 


The model presented here is based on the linear stability of a radially 
expanding fluid film. The model includes surface tension, but does not 
include viscosity. 


Geqeidae a cylinder of fluid, denoted fluid 1, of radius R,(t), where t 


represents time, and height 2z2,. This is surrounded }y an annulus of 
fluid, denoted fluid 2, of outer radius Rj (t) and also of height 2z,. 
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Surrcunding the annulus of fluid 2 is an infinite fluid medium denoted 
fluid 3. The interface between fluids 1 and 2 is denoted interface 1, 
and the interface between fluids 2 and 3 is denoted interface 2. A 
cross section of this arrangement is shown in Figure 8. 


With ref rence to a fuel-air cloud, fluid 1 represents the detonation 
reaction product gases for the high explosive burster charge. Fluid 2 
-vepresents the liquid fuel to be dispersed. Fluid 3 represents the 
ambient atmosphere. | 


Let 

represent fluid density, 

represent the radial spatial variable, | 
represent time, 

represent the radial velocity component, and 
represent pressure. 


e 
VGod wt} 


We non-dimensionalize using the following scales: 


° 2, for length, 

¢ VY. for velocity, 

¢ 2,/V, for time, and 

° (Pp, + Py + p,)V,” for pressures. 


We also define the density ratios f, = p,/(p, + Py + Pg): 
This produces the non-dimensional variables 
e r for radius, 
e 7 for time, | 
e wu for the radial velocity component, and 


¢ p for pressure. 


The correspondence between dimensional and dimensionless variables is 
given in Table 1. 


We assume that the fluids are inviscid and incompressible, and that the 
flow is purely radial, irrotational, and isothermal. 


3 IG 


Interface 2-— 


Interface 1 


Fluid 1 


Fluid 2 


Fluid 3 


Figure 8: Oylindrical geometry used for the film instability model. 
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Although these assumptions are not valid in the expanding detonation 
reaction product gases, or in the surrounding atmosphere, they permit an 
analytical solution. to be obtained. Predictions from this solution will 
provide useful bounds on the liquid breakup time and the initial drop 


diameter. These bounds can then be compared to the experimental data. 


Table 1: Relationship between dimensional variables and parameters and 
| dimensionless variables and parameters for the film instabil- 
ity model. 


Divide , | to Obtain 
Dimensional Variable ; Dimensionless Variable 


Radius, R 


Radial Velocity, U 
Time, t 
Hens ps 
Pressure, P 


Surface Tension, 0 


The basic motions of the inner and outer interfaces are prescribed by 
r,(T) and r)(T), respectively, with 


Ret) EG) £2) 


where the overdot indicates differentiation with respect to T. This 
relation is a consequence of the assumption of isochoric flow. The 
basic motions are then perturbed by small amounts 7,Y 2, to 


r; (7) + 94 (7) Yen (inner surface) 


it 


2 (7) 


ro (T) r (7) + No (T) Yen (outer surface) 


201 = 


Tn ee HI m my Ripa re Tig de apne 


where Pee A cos (m@) + BS sin(mé@) , and Z0 i cos (nt) 


A, and B, are constants. The index n specifies the axial mode; the 
index m specifies the polar (or circumferential) mode. | 


Analysis peyeals that under the assumptions: made. above, we must have the 
same values of n and m for the two interfaces, although the magnitudes | 
1 and Ne need not be equal, and generally are not. 


We find that to first order, the evolution of the disturbances 7, 
governed by the equations 


, a 
d qd {/12 4 
dr eat) - 4191 - "3 ar oe dr (0%) 0 (1) 
eee ee yg agree (ee : 
ar \"olaafe) ~ Sole * > ar |r, ar (Fa) (2) 
where 
G Bb. ~ By) v2 72 - 7 d_ = ; + F aes kr 7 
ie aa ( 2 i) ivi i dT Tr; 11 rw ( } 


Py Talks) Po Tales) Ky’ (9) = Ta’ (hte) Balhry) 
Aa“ Tey TO) OF) 1) Oy) 
fine Po T,, (kr) K,’ (kr) 7 Ty (#3) KE, (eri) 

12 RTS (ee) () ~ Ta) KO) 
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In these equations, I, is the modified Bessels function of the first 
kind of order m, and K_ is the modified Bessels function of the second 
kind of order m. A prime denotes differentiation with respect to the 
_ argument. , | | : 


We, = (p, ao Ps)V.724/0; is the Weber number of interface i. 

Equations (1) and (2) resemble a Sturm-Liouville system. However, owing 
to the coupling between ", and N,, it is not possible to assert a priori 
that the solutions are stable if the functions G, are negative, nor that 


the solutions are linearly unstable if the functions G, are positive. 
(This is because examples of such systens exist which violate this.) 


Equations (1) and (2) reduce to the appropriate result foe a single 
interface in the limit ane and to the classical result for an 
infinite plane layer in the limit ae a as 

Equations (1) and (2), with appropriate initial conditions, 7,(0) and 
N,(0), can be integrated numerically, and the behavior of the eolueions 
determined for various values of f., various disturbances (m, n), and 


various functions r,(T). Cases of interest in the explosive dispersal 
of liquids include: 


7 p. > B,; Bs, 

° Bf, > By > By; 

° constant velocity expansion, 

e decelerating expansion, and 

° accelerating-decelerating expansion. | 


The behavior of the solutions to (1) and (2) for these cases is explored 
in Appendix 4. 


We restrict our attention to (dimensional) functions R (7) of the form 
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R(t) = Raft - exp[-a(t + to)]] 


which represents a reasonable approximation, based on cloud growth data 
for typical fuel-air explosive clouds. R, is the fully expanded, 
equilibrium cloud radius obtained from experimental data, if available; 
if Ry, is unknown, 30 R, is a reasonable estimate; a is a time constant 
obtained from experimental data, if available; otherwise a value of 
0.3/msec represents reasonable estimate. t, is a constant, determined 
from | 


R (0) = R, {1 - exp[-at,)} 
In dimensionless variables, 
r(7) = ry{1 - exp[-a" (7 + T)} 


where @ = Zq/R. and ry = Ry/z,. A convenient velocity scale is then 
o 


| Examination of typical values for f, and typical expressions for ri (T) 
show that, in the absence of surface tension (see Appendix 4) 


¢ Radially expanding films are inherently unstable: 
there always exist modes which are amplified, even 
for a constant velocity expansion. 


¢ The instability may be oscillatory: the amplitudes 
behave like exp(At), where Im{A} # 0. 


° Disturbances with smaller values of n and m tend to 
be more unstable than those with larger values. 


Adding surface tension tends to stabilize the films, but does not elin- 
inate the instability. 


To use the model, we can estimate the breakup time and initial drop dia- 
meter from the time when the amplitude of the instabilities equals the 
thickness of the film, for some appropriate initial conditions »,(0) and 
9,(0), t.e., the breakup time 7, is the time when 
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where w(7) is the dimensionless film thickness, and the initial drop 
diameter is estimated to be 


ds 2w (7) : 


where d is the dimensionless drop diameter. The dimensional drop 
diameter, D(t), is d(7) Bg » and me dimensional film thickness W(t) is 
W(T) Zg. 


The linear. stability model is valid only for In,(7)| very much less than 
one. To predict the breakup time and the initial average drop dia- 
meter, it is necessary to violate this criterion; i.e., to use the 
linear model in a flow regime where a nonlinear model is required for 
accurate predictions. This means that the predicted breakup time will 
be too small, and the predicted initial arithmetic average drop diameter 
will be too large. 


‘The predicted breakup time is too small because the growth rate of these 


instabilities is known to decrease in the nonlinear regime [21]. Con- 
sequently, the estimates of the breakup time made with the linear film 
instability model will underestimate the breakup time. 


The predicted initial arithmetic average drop diameter will be too large 
because in the nonlinear regime the small-scale components of the 


evolving disturbance grow faster than the large-scale ones [21]. Thus 


the disturbance modes that cause the breakup will be smaller (i.e., will 
have larger values of n and/or m) than the critical ones indicated by 
the linear film instability model. Hence the breakup will produce 
smaller initial arithmetic average drop diameters than those indicated 
by the linear film instability node]. 


However, the model can be used to place a lower bound on the breakup 
time and an upper bound on the initial drop diameter. In the absence of 
more precise data, these bounds provide useful information for parameter 
studies with the far-field dispersal model. 


In addition, since the initial conditions are unknown, we can take the 


largest permissible values of [n,(0)| and |7,(0)| to produce the fastest 
growing disturbances (in the linear model) and hence the earliest 
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estimated breakup time and the largest estimated initial arithmetic 


average drop. diameter. 


The initial conditions used for integrating (1) and (2) were 


1, (0) = -0.01, 7,(0) = 9, (0) = 4, (0) = 0. The evolution of 9, and w(T) 


is shown in Figure 9. Numerical results are given in Table 2. The 
experimental device used by Samirant etal. was a cylindrical annulus 
with an inner radius of 1.5 cm, an outer radius of 5.5 cm, and a height 
of 12 cm, containing ethylene oxide. 


Table 2: Comparison of the breakup time and drop diameter 
: predicted by the film instability model with the 
experimental results of Samirant et al. for 
ethylene oxide. 


, Arithmetic 
Source Breakup Time | Average 
| Drop Diameter 
[msec] [cm] 


Samirant et al. [12] 


Film Instability 


* For the film instability model, the drop diameter at breakup 
was computed from d = w(7,). 


As expected, the estimated breakup time is less than the measured 
breakup time. The estimated initial arithmetic average drop diameter is 
slightly less than the measured drop diameter. However, according to 
Samirant et al. [12], immediately prior to breakup the film thickness 


was 3 mm, and immediately after breakup the arithmetic average drop: 


diameter was 5 mm. The film instability model predicts a film thickness 
and an arithmetic average drop diameter at breakup of 4.4 cm. This sug- 
gests that the linear model slightly over-estimates the film thickness 
at breakup, and that the initial arithmetic average drop diameter is 
between one and two times the film thickness just prior to breakup. 
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Figure 9: Evolution of the disturbances 7, and the film thickness 
w(T). 
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3.2 Summary of the Film Instability Model 


' ‘The film instability model places a 


¢ lower bound on the breakup time, and 


e lower bound and an upper bound on the initial 
arithmetic average drop diameter. 


To use the model, we estimate the breakup time and initial drop diameter 


from the time when the amplitude of the fastest growing instability 


ne not 


(obtained by integrating (1) and (2) numerically) equals the thickness 
of the film, for initial conditions 7,(0) and 9,(0) given by 


I7,(0)| = max{min{0.01, (1.0 mm) /Zg}, 0.01} 


and | 
|7, (0) | = max{min{O.01, (1.0 m/sec) /(aRy)}, 0.01} 


where aR, is the velocity scale. 
The dimensionless breakup time T, is the time when 
max{|1, (7,) F |73 (7) |} = r (7) ~ Ty (7,) = w(7,) 


where w(T) is the dimensionless film thickness, and the initial average 
dimensionless drop diameter d is 


w(7,) $d 2 w(7,). 
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4. THE COMBINED NEAR-FIELD DISPERSAL MODEL 


The final near-field dispersal model is a eoubanneien of the OTH hydiseus 
dynamics code and the film instability model. 


OTH provides predictions of the temperature, pressure, and velocity 
fields, and the bulk material distribution. 


Lower and upper bounds on the initial arithmetic average drop diameter 


are obtained from the film stability model, which was described in 


Section 3. 


This information is provided to the far-field model as a OTH restart 
file (at a time after the predicted breakup time), and a prediction of 
the initial arithmetic average drop diameter. 


Results from the simulation of a typical near-field dispersal are shown 
in Figures 3-7 for the initial device shown in Figure 2. The OTH input 
files and run history milestones for this simulation are given in 
Appendix 1. 


The breakup time and initial arithmetic average drop diameter at breakup 
predicted for this simulation by the film stability model are 


t= 0.45 msec and 1.4 cm § D,, § 2.8 com. 


These were calculated with 7, = 9, = -0.01, n, = 9, = 0.01, and no sur- 
face tension. Surface tension was ignored at the inner interface 
because the high temperature of the inner interface (Figure 5) indicates 
that the surface tension is very small there. The estimated Weber 
number for the outer interface is on the order of 8 x 10°; from the 
definition of the function G, (following equation 2), this large Weber 
number introduces a negligible influence on the evolution of 9,. 


The continuation of a near-field simulation into the far-field is 
described in another report [22]. 
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5.  CONOLUSIONS AND RECOMMENDATIONS 


A model has bean described for predicting the conditions in the early 
phases of the formation and expansion of a fuel-air cloud, for liquid 
fuels. In the archetypical fuel-air device, the liquid fuel is 
contained in a canister in the form of a cylindrical annulus. The 
central cylinder contains the high explosive burster charge. 


The model consists of the OTH hydrodynamics code and an ‘nalytical 
model, the film instability model. OTH provides predictions of the 
temperature, pressure, and velocity fields and the bulk material 
distribution. The film instability model provides a lower bound on the 
time at which the bulk liquid breaks up into drops and upper and lower 
bounds on the initial arithmetic average diameters of the drops. 


The model is valid in the near-field regime, where explosive forces 
dominate or are approximately equal in magnitude to the aerodynamic drag 
forces acting on the drops. Predictions of near-field cloud behavior 
and initial arithmetic average drop diameter by the model compare 
favorably with the behavior and diameter shown in high-speed movies of 
typical fuel dispersal experiments. liowever, the model remains largely 
untested against experimental data owing to the sparsity of such data. 


The model is to be used to provide initial conditions for the far-field 
model, KIVA-FAE [23], which is valid in the regime where aerodynamic 
forces on the fuel drops dominate the explosive forces. 


The near-field dispersal model ignores the effects of the canister 
fragments on the breakup process and the flow field. The fragments may 
induce turbulence in the expanding cloud; such turbulence would improve 
the comminution of the fuel and the air, producing a more uniform fuel 
distribution in the cloud. The fragments may also influence the liquid 
breakup, since at very early times the liquid squirts through the gaps 
between the fragments. 


The near-field model is inapplicable to solid fuels; an entirely 
different near-field model will be required to simulate the dispetsal of 
solid fuel particles. The degree to which solid fuels are compacted, 
agglomerate, and shatter is unknown. Simulations of explosive dispersal 
of a loosely packed powder (0.5 void fraction) and high fuel-to-burster 
mass ratios (35:1) [24] indicate that significant compaction of the 
solid particulate powders occurs. In contrast, in real fuel-air 
explosive devices the solid fuel powders will be compacted to high 
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densities when the fuel canisters are loaded, and the fuel-to-burster 
mass ratios will be higher, on the order of 100:1. The effects of these 
differences on the particle sizes and particle size distributions in the 
fuel-air cloud are unknown. 


The near-field dispersal model provides a reasonable set of initial 
conditions for the far-field dispersal model, Validation of the model 
must await the availability of experimental data in the near field 
regime of a fuel-air explosive cloud. This model is the first to 
address the behavior of expanding liquid fuel in the near-field regime, 
and in particular, it is the first model to address the breakup of 


liquid fuels into drops, 


a ay & 


APPENDIX 1: OTH Input for the Simulation of the Dispersal of Decane 
into Air 


The input files and run history for the OTH simulation of the dispersal 


of liquid decane into air by detonation of a central burster charge of 
PBX90N10 (Figures 2-7) are presented in this appendix. The files are 


presented for reference only; no attempt will be made to explain all the 


entries in these files. The interested reader is referred to the latest 
version of the OTH User’s Guide ([15]; contact the authors) for 
explanation of the structure of the files and the significance of the 
entries. 


OTH consists of a preprocessor called CTHGEN, the main computational 
module called OTH, and a file post-processor called OTHED, and a graph- 
lcs post-processor called CTHPLT. 


CTHGEN establishes a restart file for OTH. The CTHGEN input file 
specifies the computational grid, and specifies the materials and their 
location in the grid. In the cthgen input file, all characters to the 
right of an asterisk are ignored. Both CTHGEN and OTH use the ANEOS 
analytical equation of state models [15]. 


CTH performs advances the simulation in time, starting with the system 
state specified in the restart file. As in the cthgen input file, all 
characters to the right of an asterisk in the cth input file are 
ignored. | 


CTHED is an interactive program that used to examine data encoded in the 


restart file, to edit the restart file, and to create files from which. 
various one-, two-, and three-dimensional plots can be made with CTHPLT. 


Milestones in the run history for the simulation of the dispersal of 
liquid decane into air were: 


¢ The computational grid shown in Figure 2 was used from a 
problem time cf O msec to a time of 0.25 msec. 


* At 0.25 msec, the problem was rezoned to a uniform grid 
with computational cells 1 cm in the radial direction and 
1 cm in the axial direction. The small aluminum plug was 
removed during the rezoning, as it was no longer 
influencing the flow significantly. 


¢ Following the rezoning, une simulation was extended to a 
time of 2 msec. 


Total cpu time ere for the simulation was on the order of 52 cpu 
minutes. 
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cthgen input 
Saletetetae fee penn nn en pen ne pn to ee en penne nt 
Simulation 16e -- Coarse Mesh. 


This is a test FAE problem to track fuel concentration as a func- 


tion of time for a ‘typical’ NWO jug FAR device. The device is 


an idealization of an Arrowhead purified water bottle. The water 
bottle itself is ignored, as is the glass tube containing the HE 

burster charge. The fuel is decane, and the burster charge is | 
PBX9010 (90% RDX, 10% Kel-F). The Aneos input values for decane 

were devised by D. R. Gardner; the PBX9010 Aneos input values are 
the appropriate JWL parameter values. The aluminum plug atop the 
burster charge is modeled using a library equation of state. 


The PBX9010 explosive is detonated as a line charge, which extends 
from z=0 to z=12.70 cm. The charge is 0.98032 cm in radius. The 
outer radius of the decane is 13 cm. The half-height of the device 
is 17.78 cm (7 inches). The simulation is axisymmetric, and by 
virtue of symmetry, only one half of the axial length is modeled. 


Z 
| Air Air 
| aca. 
| | | 
| Al | | 
ie | 
| | 
ie | 
es | 
{3 | | 
IX | Decane { Air 
|9 | | 
|O | | 
{1 | | 
[O | | 
| | | 
| | 
_| 5 ne: 
aoe nee of ep pe en en pe ep ee pe eet 


Title 


i 
CTH Simulation 16d: Decane dispersed by PBX9010 
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Kee eee ne ce ae ee pom eeennee poe nee $m ee pee nme paneer e ene +. 


ok 


* 
. CONTROL | 
# BGTEMPERATURE = 0, 02567785 * Background beuveravixe for voids. 
MMT * Give each material its own temperature. 
* TIME = 0.0 * Problem start time. 
* VISCOSITY BL=0.1 BQ = 2. * Artificial viscosity coefficients. 
ENDCONTROL ; 
EDIT 
Block = 1 
- Noexpanded 
Endblock 
. ENDEDIT 
x, 
fe 
kee ee fee ee peer e wee ner pono perm enn freee ree + 
a 
MESH 


Block 1 Geometry = 2DCylindrical Type = Eulerian 


xO = 0.0 
x1 Number | of cells = 5 Width = 1.0 Ratio = 1.0 
x2 Width = 19.0 DXFirst = 0.2 DXLast = 1.0 
x3 Number_of cells = 80 Width = 80.0 Ratio = 1.0 
Endx 
‘ | 
yO = 0.0 | | 
yl Number of rare: 100 Width = 100.0 Ratio = 1.0 
Endy 
x 
* Activate small portion of mesh. 
* 
Xaction = 0.0, 5.0 
Yaction = 0.0, 15.0 
* 
Endblock 
ENDMESH 
eo 
Ko ee paren ene +---------+ tome nee tome meee pee e eee pone en nee + 


Insertion_of Material 
* : 
Block 1 
- 
* 


535° 


PBX9010 


Package = 
Material = 3 
Density = 1.787 
Numsub = 10 
Pressure = 1.01325E+6 © 
* \ n 
: Insert Box : | 
xl = 0.00 x2= 0.98032 . 
yl = 0.00 y2 = 12.70 
| Endinsert _ 
* 
-. Endpackage 

* 7 

* : : 

Package = ’Aluminum Plug’ 
Material = 4 
Pressure = 1.01325E+6 
P : 
Insert Box _ . 
xl = 0.00 x2 = 0.98032 
yl = 0.00. y2 = 17.78 
Endinsert 
* 
' Endpackage 

‘ | 

7 | | 
Package = ’Decane’ 
Material = 2 
Density = 0.7263 
Pressure = 1.01325E+6 

* 

Insert Box | 

xl = 0.00 x2 = 13.00 
yl = 0.00 y2 = 17.78 
Endinsert 

* | 
Endpackage 

* 

: | 
Package = ’Dry Air’ 
Material = 1 
Pressure = 1.01325E+6 

* 

Insert Box 
x1 = 0.00 x2 = 100.0 
yl = 0.00 y2 = 100.0 
Endinsert 

R | 
Endpackage 

* 

* 


an | 1 rams q ' ‘ow " ‘ » 
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* 


Endblock 


_Endinsertion ~ 


* 
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* 


* 
* 


* Equation of State Information. 


cy ae 


* my 
EOS Number_of materials = 4 
: = 
Materiall = -1 * Dry Air 
Material2 = -2 +* Decane 
Material3 = -3 * PBX9010 
Material4 = -4 * Aluminum 
* 
iv 
aneosO -1 ’Dry Air’ | ; 
aneosl 3.0 -2.0 1.1830E-3) 0.0 0.0 0.0 0.400 0.000 
aneos2 O 0.0 8.33924E10 0.0 0.0 00 O 0 
aneos3 0.0 0.0 0.000 0.0 0.0 0O 0 0 
aneos4. 7.0 0.7809 | 
aneosS 8.0 0.2095 
aneos6 18.0 0.0096 | 
* : : 
*ANEOS inputs for decane were compiled by D. R. Gardner (6425). 
aneosO -2 ’Decane’ | 
aneosl 2.0 3.0 0.7700 0.02097742 1.01325E6 1.500E10 0.300 
0.000 | 
aneos2 O 2.0 9.42E10 0.02097742 0.0000 0.0 O 0 
aneos3 0.5 ©@.0 0.000 0.00 0.000 0 0 0 
aneos4 1.0 22.0 
aneoso 6.0 10.0 
rs | , 
*ANEOS inputs for pbx9010 were compiled by D. R. Gardner (6425). 
'-aneosO -3 ’PBX9010’ | | 
aneosl 6.0 -3.0 1.7870 0.02567785 1.01325E6 1.0 0.350 0.0 
-aneos2 0.0 0.0 0.000 0.0411 0.0000 0.0 0O 0 
aneos3 0.0 5.814 0.06801 4.10 1.00. 0 0 0 
aneos4 1 0.2694 
aneosS 6 0.1537 
aneos6 7 0.2694 
aneos7 8 0.2694 
aneos8 9. 0.0258 
aneos9 17 0.0095 | 
* : 
aneosO -4 ’Aluminum’ LIB = 3 TYPE = 0 
* 
ENDEUS 


 HEBURN 


* 
Material = 3 | 
Detonation Velocity = 8.500E+5 — 
Rho = 1.787 = 
Temperature = 0.02567785 
Number_of_cells_in_burn_front = 2 
Predetonation Pressure = 5.0B+15 


RP RP RP RP RP 


‘Dline = 0.0, 0.0 to 0.0, 12.70 
Radius = 0.98032 
Time = 0.0 : 


RP BP 


Ok 
-ENDHEBURN 


« 


TRACER its , 

* Add Lagrangian tracer particles. 

* ; 

* On each side of the PBX9010-decane interface: 
ADD 0.90, 2.0 TO 0.90, 10.0 NUMBER 5 
ADD 1.10, 2.0 Tu 1.10, 10.0 NUMBER 5 

* : ; 

x On the decane side of the decane-air interface: 
ADD 12.5, 2.0 TO 12.50, 10.0 NUMBER 5 


* 

ENDTRACER | | 

* : 

ENDINPUT * End of cthgen input 


% 
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* | ' eth input 

i. : . ; 
kee pone nen pone ee ee fone ee ee fee mn + 
* ‘ f 

x Simulation 16d -- Coarse Mesh. 

* | | 

x This is a test FAE problem to track fuel concentration as a func- 
* tion of time for a ‘typical’ NWC jug FAE device. The device is 

* an idealization of an Arrowhead purified water bottle. The water 
x bottle itself is ignored, as is the glass tube containing the HE ~ 
* burster charge. The fuel is decane, and the burster charge is 

+ PBX9010 (90% RDX, 10% Kel-F). The Aneos input values for decane 

* were devised by D. R. Gardner; the PBX9010 Aneos input values are 
+ the appropriate JWL parameter values. The aluminum plug atop the 
* burster charge is modeled using a library equation of state. 

| 

* The PBX9010 explosive is detonated as a line charge, which extends 
* from z=0 to z=12.70 cm. The charge is 0.98032 cm in radius. The 
* outer radius of the decane is 13 cm. The half-height of the device 
* is 17.78 cm (7 inches). The simulation is axisymmetric, and by 

* virtue of symmetry, only one half of the axial length is modeled. 
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_ CTH Simulation 16d: Decane dispersed by PBX9010 
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APPENDIX 2: A Surface Tension Model 


Oonsider a cylindrical annulus of fluid. The fluid is assumed to be 
incompressible. The length of the annulus is 224; its inner radius is 
R,(t) and its outer radius is Ro(t). The surface tension at the inner 
and outer surfaces is 9. We assume that the fluid is moving radially 
only, and that the flow is isothermal. The free energy of the surfaces 
is 


Eo = 2nz40 [R, (t) + R(t)). 


Suppose that the film breaks up into a monodisperse field of drops with 
diameter d. If the total volume of the liquid is V, then the number of 
drops N and the total surface area of the drops A, is 


‘N= 6V/(rd®), and A, = Nyrd? = 6V/d, 
The total free energy of the drop surfaces is 


E, = oA, = 60V/d. 


The total volume of fluid is V = 2z_9m [R, 7 (0) - R,7(0)]. 


Initially, and for some 0 < t < t,, the free energy of the annular 
surfaces is less than the free energy of that same fluid if it existed 
as drops of diameter d: E, < E,. As the fluid expands, it reaches a 
state where the two free energies are equal: EL = E,. At later times 
the free energy of the annular surfaces is greater than the free energy 
of that same liquid if i+ existed as drops of diameter d: E, > B,. At 
these later times, it is energetically more favorable for the fluid to 
exist as drops, rather than as a film. We assume then that breakup of 
the film into drops occurs when the two free energies are equal: E, = 
E,- 
We can thus define a parameter Ewa = E,/E, such that breakup occurs when 
E_, = 1. This leads to the expression: : 


E, fd(r.(r) + 3, (7)] | 
= = = 1 (2.1) 
E. 1 - (R,/R,)* 


where r,(T) = R, (t) /Z,, ro(7) = Ri (t)/2,, 7 = t/T,, RB, = R,(0), Ro = 
R.(0), d= d/z,, and f = (z,/R,)*/3. T, is a time scale. 
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Equation (2.1) can be solved for d: 


i t= RR 
Aiea: eee as Ss (2.2) 
po ri(7) + ry(7) 


Equation (2.2) can be evaluated if expressions are known for r,(T) and 
oe 


Assuming that the fluid is incompressible and moves only in the radial 
direction allows us to write 


Qnz_(R,7(t) - R,"(t)] = 2nz,(R,"(0) - R,”(0)] 


and hence Rj (t)U,(t) = Ry(t)U,(t), where U is the velocity. Then 
assuming that R(t) has the form 


R(t) = Ref 1 - exp[-a(t + t,)]} (2.3) 


we can derive expressions for r.(T) and r,(7). Equation (2.3) 
represents a reasonable approximation, based on cloud growth data for 
typical fuel-air explosive clouds. R, is the fully-expanded, 
equilibrium cloud radius, @ is a time constant, and t, is determined 
from 


R,(O) = Ry{ 1 - exp[-at,]} 


From (2.3) we can construct a time scale using aR, and 2,: T 
Z4/ (aR). Then 


r,(T) = Te{ 1 - exp[-a"(7 + T,)]} (2.4) 


where @, = Z3/Ry and ry = Roo/Zy. 


Equations (2.2) and (2.4) provide an estimate of the initial drop 
diameter at breakup, if the breakup time is known. Conversely, if the 
initial drop diameter is known, then the breakup time can be estimated. 


Predictions from the model can be compared to the measured breakup time 
and initial drop diameter reported by Samirant et al. [12] for a device 
in which 2, = 6.0 cn, R\(O) = 5.5 cm, and R,(0) = 1.5 cm. Using the 
velocity data presented in [12], R, = 184.10 cm, t, = 0.020641 msec, and 
a = 1.4695/msec; thus a” = 0.032591, r, = 30.683, T, = 0.93069, and 


e ASS 


r, (7) = 30.683 { 1 ~ exp[-0,032591(r + 0.93069) }} 


Predictions for the initial arithmetic average drop diameter as a 
function of breakup time for this model device are given in Figure 10. 


Predictions of the breakup time as a function of the initial arithmetic 
average drop diameter for this model device are given in Figure 11. 


Predictions of the initial arithmetic average drop diameter and breakup 
time are given in Table 3, and compared with the Meee: values of 
Samirant et al. [12]. 


The predicted breakup time (when the initial drop diameter is known) is 
too small by a factor of three. The predicted initial drop diameter 
(when the liquid breakup time is known) is too small by a factor of two. 
The predicted values are of the correct order of magnitude, but do not 
provide very accurate predictions of the liquid breakup time and the 
initial drop diameter. However, considering the assumptions made in the 
model, this agreement is quite good. 


Since either the breakup time or the initial arithmetic average drop 
diameter at breakup must be known to use this model, and since in 
general neither of these quantities is known, this model is not 
especially useful in dispersal modeling for fuel-air explosives. 


Table 3: The surface tension model predicts breakup 
quantities that compare favorably with the 
experimental results. 


Arithmetic 
Average 
Source Breakup Time | Drop Diameter 


Samirant et al. [12] 


Surface tension 
(assuming drop dia- 
meter of 0.5 cm) 


Surface tension 
(assuming breakup 
time of 1.3 msec) 
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Figure 10; The initial arithmetic average drop diameter predicted by the 
surface tension model, as a function of the breakup time. 
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Figure 11: The breakup time predicted by the surface tension model, as a 
function of the initial arithmetic average drop diameter, 


APPENDIX 8: ‘The Grady Fracture Model 


_ The Grady fracture model was developed by D. EB. Grady (1534) for esti- 


mating fragment sizes in the explosive fracture of solids, It has been 
used successfully to estimate the sizes of fragments produced in the 
brittle fracture of oil shale and steel [16], the explosive loading of 
cylindrical shells of uranium [17], the impact fragmentation of lead and 
uranium plates [18], and the brittle and ductile spall of a variety of 
metals [19]. 


The assumptions underlying the Grady model do not limit it a priori to 


solids. However, comparison with the limited data available for liquid 
spall suggest that is not as useful as the models presented in Section 3 
for estimating the breakup time and initial drop diameter for 
explosively dispersed liquids. 


The Grady model is described here briefly, ind some predictions are made 
with it for comparison with the data of Samirant et al. [12]. The 
predictions of the Grady model do not match this data as well as the 
predictions given in Section 3. 


The Grady Fracture Model 


As applied to liquids [19], the Grady fracture model takes the following 
form for energy-limited spall (see also Table 4): 


For spall dominated by surface tension: 


1/2 1/3 | 1/3 
as 1 (60 480 
P, = |2p 0 we j t, = | d = | . 
b ( O bo) pee b pe" 


For spall dominated by viscous dissipation: 


evi/2 . 1/2 *1/2 
Le [2p 0° ye] » ty» [2u/pce é] » and D, = (8u/pe) / 


In these equations, P, is the spall strength of the liquid, and t, is 
the time at which the liquid breaks into ‘‘fragments’’ of size D. o is 
the liquid surface tension, p is the liquid density, 0 is the speed of 
sound in the liquid, wis the dynamic viscosity of the liquid, and é is 
the strain rate. 


me yee 


Using the previous approximation for the location of the outer boundary 
of the expanding liquid shell, we have 


€ = (vB/a){exp[-a(t+t.)] - 1}, 


The breakup time and initial arithmetic average drop diameter predicted 
by the Grady fracture model for the device used by Samirant et al. [12], 
using the values of the material properties for ethylene oxide shown in 
Table 5 (in the absence of a known value, co has been approximated using 
a value typical for similar liquids) are given in Table 6. Included are . 
predictions for spall dominated by surface tension or by viscous 
dissipation. | 


For either extreme, the breakup time and initial arithmetic average drop 
diameter predicted by the Grady model are both far too small, compared 
to the experimental data. This is probably because in the explosive 
dispersal of liquids, spalling is not the dominant mode of liquid 
breakup [1, 20]. 


Table 4: The equations for the Grady fracture model. 


‘Parameter Surface Tension Viscous Dissipation 


Dominated Spal] 


Dominated Spall 
Breakup Time 
t 
b 


Drop Diameter, 
Dp 


Spall Strength, 
Py 
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Table 5: ‘Ethylene oxide property values. 


Density, kg/m 


Surface Tension, N/m 2.42 x 10° 


Dynamic Viscosity, kg/m-s 


Sound Speed, m/s | ~ 1.6 x 10° 
(estimated from similar liquids) 


Table 6: Breakup time and initial arithmetic average drop diameter 
predicted by the Grady Fracture Model. 


Parameter Breakup Time, t, | Initial Average 
[msec] Drop Diameter, d, 


[em] 


Surface Tension i 
Dominated 2.1x 10 


Viscous Dissipation 
Dominated 


Data from Samirant 
et al. [12] | 


AG = 


APPENDIX 4: Development of the Film Instability Model 


Presented in this appendix are some details of the development of the 


film instability model (Section 3), and the parametric study of the 


model for cases of interest in fuel-air dispersal studies. 


Consider a cylinder of fluid, denoted fluid 1, of radius R, (t), where t 

represents time, and height 2z,. This is surrounded by an annulus of 

fluid, denoted fluid 2, of outer radius R,(t) and also of height 2z,. 

Surrounding the annulus of fluid 2 is an infinite fluid medium denoted 

fluid 3. The interface between fluids 1 and 2 is denoted interface 1, 

and the interface between fluids 2 and 3 is denoted interface 2. A 
cross section of this arrangement is shown in Figure 8. 


With reference to a fuel-air cloud, fluid 1 represents the detonation 
reaction product gases for the high explosive burster charge. Fluid 2 
represents the liquid fuel to be dispersed. Fluid 3 represents the 
ambient atmosphere. 


Let 

© p represent fluid density, 

° R represent the radial spatial variable, 
represent time, 
e U represent the radial velocity component, and 
e P represent pressure. — 


e 
ne a 


We non-dimensionalize using the following scales: | 


° gz, for length, 

¢ V. for velocity, 

° 2z,/V, for time, and 

© (Pp, + Pe + pov for pressures and pressure 
differences. 


We also define the density ratios f, = p,/(p, + Py + Pg,)- 
This produces the non-dimensional variables: 

e r for radius 

e 7 for time 


e u for the radial velocity component, and 
e p for pressure 
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We assume that the fluid in the film is inviscid and incompressible, and 
the flow is purely radial, isochoric, irrotational, and isothermal. 


Equivalently, we can assume that the flow is purely radial, the 
Mach number M + 0, the Reynolds number R, + ©, and the Prandtl 
number P_ + 0 (but. in such a way that RP. + 0), 


Although these assumptions are not valid in the expanding detonation 
reaction product gases, or in the surrounding atmosphere, they permit an 
analytical solution to be obtained. Predictions from this solu.ion will 


“provide useful bounds on the liquid breakup time and the initial drop 


diameter. The predicted bounds can then be compared to the experimental 
data. | 


The basic motions of the inner and outer interfaces are prescribed by 
r,(T) and r,(7T), respectively, with 


(1) $4 (7) = rE(r) F(7) 


where the overdot indicates differentiation with respect to T. This 
relation is a consequence of the assumption of isochoric flow. 


The basic motions are then perturbed by small amounts 7,Y 2, to 


ry (7) r, (7) + 1, (7) Yen | (inner interface) 


To) (T) ro (7) + No (T) Yen (outer interface) 


: mea Yo = A cos (m6) + B Sin (md) , and Z = cos(nWz). 


A. and B| are constants. The index n specifies the axial mode; the 
index m specifies the polar (or circumferential) mode. 


We have assumed here that the indices (m, n) are the same for both 
interfaces. It can be shown by repeating this entire analysis with 


different indices for each interface that the only solutions with 
different values of the indices for each interface are trivial ones. 


on -) a 


an] rar py 1 oo | ! my nou A " ' hort oa f yore hone Trae a yy 


Herpes ert 


Hence it is valid tor assume that the iad twee. are the same foe both 
interfaces. 


Under these conditions, the equations of motion in the three fluids 
admit a solution in terms of a velectey potential ¢, in each ead such 


that 


+ 
Voges = O and v= Ve. 


Vis the vector gradient ee and V* is the Laplacian. 


As boundary conditions, we require that the normal velocities of the 
fluids be equal at the appropriate interfaces: 


+ F a e a a 
n, is a unit normal to interface i. 


We follow the method used by Plesset [25]. In Fluid 1, we find that, to 
first order in 9,, $, is given by : 


Abts Mg te@)|: Ter) 
4 7 r; (7) r,(T) In r+ 1 sc (r)| k 1. *(kr.) 


O<r <7, (7), and k = 


I_ is the modified Bessels function of the first kind of order m. A 
prime denotes differentiation with respect to the argument (kr in this 
case). 

In Fluid 3, we find that, to first order in 7,, ¢, is given by 


| : foo. fe) 2. K (ke) 
#3 - r (7) r (7) In r + No CE) KK (ke) ‘n4n ’ 
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| ro (7) <ece 


| K is the modified Bessels function of the second kind of order n. 


In the film, we find that, to first order in 7,, $, is given by 


to = 7, (7) Pada r+ myo 


es ee es 
1 73 1 To n 


I (kr) K (kr) - I,’ (kr). K,, (kr) 
a a ar rT ’ 
: kT, (kr) K (kr) - I, (kr,) KO (kr) ] 


2° ETE (er) Key) ey) KO] 


Using Bernoulli’s principle, 


ee ees 
or * 2 ("#i) * 5, =e) 


in each fluid, and the interface conditions 


s ~ ot Jot. oh " _ fi 
Pus = Py OW Ps 2 (YnAn : Yn Zi ven } 


ee ee eee ea ee ak 
Pos ~ Po.” W r 2 ( mn m on mn } 
e2 | o Ty | : 
where We, = (p, + Po + Ps)V, 04/0; is the Weber number of interface i, 


and P.. ie the pressure just inside (-) interface i or just outside (+) 
interface i, we find that, to zeroth order in Na? 


aa: 


" lees 


. if d f_ 1 92 7 . 43 
(Po - A3){1n(r3 a (ital +3 fi soe Aya oye 


el 
(4:1) 
d fs 1 92 | | 1 
(Ps Pa)\ 0(*, a Fetal? 2 70 = Pgeg(t). ~ Pata(t) + FW 


(4.2) 


‘These equations relate r, and r, to the functions c, i(T). We will assume 


that, given functions r (7) and ro(7), we can find functions c, (7) 
consistent with (4,1) and (4.2). 


To first ee in ,, we Rave 


: qd {fio q : 
dr GHA ~ S19) - Ty ar rar 72%) ae ee 
d 7 a [fai | 
dr [rfoate) - Sate + To SF r,t nl me (4.4) 


where 


Paes < ges 


i 
ee een 


Po 1, {kr } K(k) - 1 (kr,) K(k) 
fon TSC) RO) — Fa Fe) Be id 
op bs K (kr) Po I, (kr) Ky (kr;) _ Th’ (Ey) K (kr) 
‘09 | : k Ko (kr) ° k I (kr) Ko (kr) 7 I (kr; ) KD (kr) 


Equations (4.3) and (4.4) resemble a Sturm-Liouville system. However, 
owing to the coupling between 7, and 7,, it is not possible to assert a 
priori that the solutions are stable if the functions G, are negative, 
nor that the solutions are linearly unstable if. the functions G. are 


positive. (This is because examples of such systems exist which violate 


this.) 


Equations (4.3) and (4.4) reduce to the appropriate result for a single 
interface in the limit r, + r,, and to the classical result for an 
infinite plane layer in the limit r, +, r+ mr, <r [26]. 


Equations (4.3) and (4.4), with appropriate initial conditions, "; (0) 
and 7’,(0), can be integrated numerically, and the behavior of the 
solutions determined for various values of f,, various disturbances (m, 
n), and various functions r,(7). Cases of interest in the explosive 
dispersal of liquids include: 


* Py > Pyr Bs: 

° By > Ba > Bg, 

* constant velocity expansion, 
¢ decelerating expansion, and 


* accelerating-decelerating expansion. 


We restrict our attention to (dimensional) functions R(T) of the form 
Ri(t) = Ry{ 1 - exp[-a(t + t,)]} 


which represents a reasonable approximation, based on cloud growth data 
for typical fuel-air explosive clouds. R, is the fully-expanded, 
equilibrium cloud radius, a@ is a time constant, and t, is a constant, 
determined from 
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R,(0) = Reo{ 1 - exp[-at,]}. 
In dimensionless variables, | | 
| ro(T) = te{ 1 - exp[-a"(r + Te) 13 


od ee = Zy/Rg and rey = Ry/z,. A convenient velocity scale is then 
= a , 
re) 0° 


To use the model, we can predict the breakup time and initial drop dia- 
meter from the time when the amplitude of the instabilities equals the 
thickness of the film, for some appropriate initial conditions 7,(0) and 
7,(0), i.e., the breakup time 7, is the time when 


max{ iy (75) | » |Yal(T») |} = Pot) ~ Fale) = CH) 
and the initial drop diameter is predicted to be w(7,). 


The model is formally valid only when |,| is very much less than one. 
From the construction of the physical device, the amplitudes of initial 
disturbances are expected to be small, on the order of a millimeter or 
so. Thus a practical limit is |7,(0)|z, ¢ 1 mm. So for computational 
purposes we take 


I7,(0)| $ max{min{0.0]1, (1.0 mm)/z,}, 0.01}. 


Similarly, there are bounds on the velocities 9,(0). A practical limit 
is disturbance velocities on the order of one meter per second; thus 


17, (0) | $ max{min{0.01, (1.0 m/sec)/(aR,)}, 0.01} 
where aR, is the velocity scale. 


As noted above, the linear stability model is valid only for |, (7) | 
very much less than one. To predict the breakup time and the initial 
arithmetic average drop diameter, it is necessary to violate this 
criterion; i.e., to use the linear model in a flow regime where a 
nonlinear mode] is required. This means that the predicted breakup time 
will be too small, and the predicted initial arithmetic average drop 
diameter will be too large. 
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The predicted breakup time is too small because the growth rate of these 
instabilities is known to decrease in the nonlinear regime [21]. Oon- 
sequently, the estimates of the breakup ae made with the linear model 
will underestimate the breakup time. 


The predicted initial arithmetic average drop diameter will be too large 
because in the nonlinear regime the small-scale componerts of the 
evolving disturbance grow faster than the large-scale ones [21]. Thus 
the disturbance modes that cause the breakup will be smaller (i.e., will 
have larger values of n and/or m) than the critical ones indicated by 
the linear model. Hence the breakup will produce smaller initial 
arithmetic average drop diameters than those indicated by the linear 
wodel , 


However, the model can be used to place a useful lower bound on the 
breakup time and an upper bound on the initial drop diameter. 


In addition, since the initial conditions are unknown, we can take the 
largest permissible values of |n, (0) | and |, (0) |. to produce the fastest 
growing disturbances (in the linear model) and hence the earliest 
estimated breakup time and the largest estimated initial arithmetic 
average drop diameter. 


Predictions from the model can be compared to the measured breakup time 
and initial drop diameter reported by Samirant et al. [12] for a device 
in which 2, = 6.0 cm, Rj(O) = 5.5 cm, and R,(0) = 1.5 cm. Using the 
velocity data presented in [12], R, = 184.10 cm, t, = 0.020641 msec, and 
a = 1,4695/msec; thus a” = 0.032501, r,_ = 30.683, T, = 0.93069, and 


r,(T) = 30.683 {1 - exp[-0.032591(r + 0.93089) ]} 


The initial conditions used for integrating (4.3) and (4.4) were 7, (0) = 
-0.01, 7,(0) = 9, (0) = 9,(0) = 0. The evolution of 9, and w(T) is shown 
in Figure 9. Numerical results are given in Table 2. 


As expected, the estimated breakup time is less than the measured 
breakup time. The estimated initial arithmetic average drop diameter is 
slightly less than the measured drop diameter. However, according to 
Samirant et al. [12], immediately prior to breakup the film thickness 
was 3 mm, and immediately after breakup the arithmetic average drop 
diameter was 5 mm. This suggests that the linear model slightly over- 


ane: aes 


estimates the film thickness at breakup, and that the initial arithmetic 
average drop diameter is approximately twice the film thickness just 
prior to breakup. | 
The conclusion is that the film instability model places a 

* lower bound on the breakup time, and 

e lower bound on the initial arithmetic average drop 


diameter. 


We can obtain an upper bound on the initial arithmetic average drop 
diameter by taking | 


d $ 2w(r,). 


where w(T) is the dimensionless film thickness, and d is the dimension- 
less drop diameter, | 


To use the model, we estimate the breakup time and initial drop diameter 
from the time when the arvlitude of the fastest growing instability 
(obtained by integrating (4.3) and (4.4) numerically) equals the 
thickness of the film, for initial conditions 9,(0) and #,(0) given by 


In, (0) | = max{min{0,01, (1.0 mm)/z,}, 0.01} 
and 


[7,(0) | = max{min{0.01, (1.0 m/sec) /(aR,»)}, 0.01}, 


where aR, is the velocity scale. 


The breakup time T, is the time when 


max { In (7) | Ing (7) |} = r (7) 7 r; (7) = w(7)), 
and the initial arithmetic average drop diameter d lies in the range 


w(7,) $d $ 2w(7,). 
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Parameter Study of the Film Instability Model 


The behavior of the film instability model was explored for a variety of 
cases of interest in the modeling of fuel-sir explosives. (ases 
considered were | | 


¢ various values of the indices n and m, 

¢ various relative values of the density ratios f,, 
¢ various forms for the basic motion r,(7), and 

* various values for the surface tensions 0,. 


The specific values for the results reported here are given in Table 7, 


The definition of stability here is less restrictive than the usual 
mathematical definition of stability. Since the differential equations 
(4.3) and (4.4) are non-autonomous (i.e., the coefficients contain the 
time variable 7 explicitly, via r, and r,) the mathematical stability of 
the equations cannot be rigorously inferred from the eigenvalues. What 
is required in the film instability model is that the amplitude of some 
disturbance mode with indices (n, m) grow until it equals the thickness 
of the film. For example, an oscillating disturbance with constant 
amplitude, which is mathematically stable, will cause the breakup of the 
film when the film thins to the point where its thickness equals the 
amplitude of the oscillations. (This means that any non-damped 
disturbance will ultimately disrupt the film, as the film thins.) 


In the film instability model, a radially expanding film is said to be 
(linearly) unstable if there exists one disturbance mode whose amplitude 
at some time equals or is greater than the thickness of the film at that 
time. 
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It would be convenient to infer the mathematical stability of a given 
disturbance mode from the eigenvalues of (4.3) and (4.4). Although this 
cannot be done rigorously, the eigenvalues will be used to indicate the 
stability of the expanding film as usual: if any one of the eigenvalues 
has a positive real part, then the film is indicated as unstable by the 
eigenvalues, This is valid for two reasons} | 


* in all cases tested, if the film was unstable (as 
defined above), then there was at least one eigenvalue 
with a positive real part, and 


° the conclusions drawn from the stability of the film 
are the same as those drawn from the stability of a 
radially expanding single interface, for which the 
eigenvalues do rigorously indicate the (mathematical) 
stability. 


In using the film instability model, the primary interest lies in when 
the amplitude of a disturbance equals or exceeds the thickness of the 
film. Thus the quantitative predictions made by the model are 
unaffected by the definition of stability. However, in order to 
understand the general behavior of a radially expanding fluid film, it 
is convenient to examine the eigenvalues. In this context, the critical 
eigenvalue is the one with largest real part, 


Plots of the largest eigenvalue for representative cases from Table 7 
are given in Figures 12-19. 


Figures 12-13 show the real part of the critical eigenvalue for various 
combinations of the disturbance indices (n, m) and for a constant 
velocity radial expansion, and for two different combinations of the 
density ratios f, of interest in the dispersal of liquid fuels. Since 
the real part of the critical eigenvalue is always positive in these 
plots, we can infer that, even when there is no acceleration, radially 
expanding fluid films are unstable. 


Figures 14-15 show the real part of the critical eigenvalue for various 
combinations of the disturbance indices (n, m) and for a constant-accel- 
eration radial expansion, for two different combinations of the density 
ratios f,. Since the real part of the critical eigenvalue is always 
ultimately positive in these plots, we can infer that accelerating, 
radially expanding films are unstable. A similar result holds for 
decelerating films (Figures 16-17). 
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Figures 18-19 show the real part of the critical eigenvalue for various 
combinations of the disturbance indices (n, m) and for an accelerating- 
density ratios f,. Since the real part of the critical eigenvalue is 
always positive in these plots, we can infer that accelerating- 
decelerating, radially expanding films are unstable. 


From Figures 12-19, we can also infer that, as a general rule, the (n, 
m) = (1, 0) disturbance is the most dangerous, i.e., grows most rapidly, 
because in most of the cases shown, its critical ‘eigenvalue has the 
largest real part. 


Although not shown in the figures, in some cases the critical eigenvalue 


is a complex conjugate pair, indicating an oscillatory disturbance with 
growing amplitude, 


The effect of adding surface tension is to decrease the amplitudes of 
the disturbances 7,. This delays the film breakup, but does not prevent 
it: the films are always ultimately unstable. The concomitant 
predicted arithmetric mean drop diameter is thus decreased. This is 
illustrated by the predicted film breakup times and concomitant 
predicted drop diameters shown in Table 8. The results in the table are 
for 


ri(T) = re{l - exp[-a"(T + 7,))} 


g = 0.97589, w(0) = 0.5 (these values give the 
same initial conditions as those for the r,(T) functions listed in Table 
7); B, = = fs = 0, p, = 1; (n,m) = (1,0); and 7, (0) = 0.01, 1, (0) = 0.01, 

7,(0) = 7,(0) = 0. The Weber numbers were W., = © (corresponding to 
parc surface tension) and various values of W.,. Decreasing W.,, and W,, 
(which corresponds to increasing surface tension) for this form of r,(T) 
produce relatively small decreases in the critical eigenvalue; the real 
part of the critical eigenvalue remains positive, indicating that the 
films are unstable. Changes in the functions G, with decreases in the 
Weber numbers are more pronounced, and are manifested in the behavior of 
Ng(7) : This behavior is illustrated in Figure 20. For W., =, n, 
increases monotonically and rapidly; while for W., = 500, 7, increases 
more slowly and oscillates. 


with r, = 21, a” = 0.05, 7 


In summary, examination of typical values for f, and typical expressions 
for r,(T) indicate that, in the absence of surface tension 
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e Radially expanding films are inherently unstable: 
there always exist modes which are amplified, even 
for a constant velocity eypansion. 


¢ The instability aay be oscillatory: the amplitudes 
behave like exp(At), where Im{\} may be non-zero. | 


¢ Disturbances with smaller values of n and m tend to 
be more unstable than those with larger values. 


Adding surface tension tends to stabilize the films, but does not elim- — 
inate the instability. | 
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Table 8: Effect of surface tension on the breakup time and initial - 
arithmetic average on diameter predicted by the film 
instability model We od, ids d.) 


Outer Interface Dimensionless Dimensionless 
Weber Number Breakup Time : Drop Diameters 
Lot 


0.00414 | 0.00828 
0.00820 | 0.0164 


0.00854 {| 0.0171 


0.00969 | 0.0174 


0.00905 0.0181 
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Figure 12: 
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The real part of the most critical eigenvalue | for four 
different disturbances for a constant velocity radial 
expansion r,(7) = 7 + 1, with density ratios A, = 0, fp, = 1, 
and f, = 0. | 
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Figure 13: The real part of the most critical eigenvalue A. for four 
different disturbances for a constant velocity radial 
expansion r)(T) = 7 + 1, with density ratios f, = 0.66667, 
A, = 0.33333, and f, = 0. | 
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. Figure 14: The real part of the most critical eigenvalue \. for four 
different disturbances for a constant acceleration radial 
expansion r(7) = 0.02577 + 7 + 1, with density ratios p, = 
0, 6, = 1, and f, = 0. 
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Figure 15: 
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The real part of the most critical eigenvalue \ for four 
different disturbances for a constant acceleration radial 
expansion 2) = 0.02577 + T + 1, with density ratios By = 
0.66667, b. = 0.33333, and By = 0: | 
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Figure 16: The real part of the most critical eigenvalue A, for four 
different disturbances for a constant deceleration radial 
expansion ro(T) = a 0257" + T + 1, with density ratios B, = 


0, f, = 1, and f, = 
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The real part of the most critical eigenvalue dh, for four 
different disturbances for a constant deceleration radial 
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The real part of the most critical eigenvalue \. for four 
different disturbances for a acceleration-deceleration 
radial expansion r,(T) = -0.016667r° + 0.02577 + 7 + 1, with 
density ratios J, = 0, 6, = 1, and A, = 0. 
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Figure 19: 
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The real part of the most critical eigenvalue 4, for four 
different disturbances for a acceleration- deceleration 
radial expansion r.(T) = -0. 0166677" + 0.0257? + 7 + 1, with 


density ratios f, = 0.66667, f, = 0.33333, and f, = 0 
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, Figure 20: Evolution of the disturbance 9, for different Weber numbers 
| Wag: Initial conditions are given in the text. w(t) is 
shown for reference. 
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